Model for erosion-deposition patterns 
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We investigate through computational simulations with a pore network model the formation of 
patterns caused by erosion-deposition mechanisms. In this model, the geometry of the pore space 
changes dynamically as a consequence of the coupling between the fluid flow and the movement of 
particles due to local drag forces. Our results for this irreversible process show that the model is 
capable to reproduce typical natural patterns caused by well known erosion processes. Moreover, we 
observe that, within a certain range of porosity values, the grains form clusters that are tilted with 
respect to the horizontal with a characteristic angle. We compare our results to recent experiments 
for granular material in flowing water and show that they present a satisfactory agreement. 



I. INTRODUCTION 

Nature preservation and environmental protection are 
important issues on the agendas of governments and non- 
governmental organizations. Consequently, these issues 
have been the subject of intense study, where the main 
concern is to understand how the actions of human be- 
ings affect the environment. For example, deforestation 
and pollutant emissions are related to climate changes 
resulting in floods and erosion. In particular, erosion can 
be responsible for diminishing the quality of life, because 
it affects the soil causing a negative impact on the econ- 
omy. Aside from its economical and ecological aspects, 
the erosion problem also attracts the interest of geolo- 
gists and physicists. In geology, this is an extremely rich 
area as many of the patterns observed in nature stem 
from erosion or deposition processes. In physics, the for- 
mation of such patterns span a huge range of spatial and 
temporal scales. This pattern formation process is di- 
rectly related to the transport of solid granular particles 
via a fluid and presents a richphenomenology along with 
a variety of applications P, Q . Particular applications 
are fractal river basins, meandering rivers, dune fields, 
granular avalanches, and ripple marks on sand banks or 
on coastal continental platforms. 

It is notoriously difhcult to provide a fully consistent 
description of particle laden flows either from a one-phase 
or a two-phase point of view. Most of the practical knowl- 
edge of erosion comes from empirical laws often derived 
from field measurements. This has provoked interest in 
theoretical descriptions of these systems [1, |3j [E| and vi- 
sualization of them in computational simulations. Sev- 
eral attempts to understand the dynamics of river basin 
formation from the statistical physics point of view have 
been made recently, but many questions are raised when 
one tries to relate basic transport properties to large-scale 
pattern forming instabilities A fundamental open 
question is the following: how do objects made of gran- 
ular materials respond to the action of external factors? 
Many experiments 0, i, i, El HH, [H, [H, [11 have been 
performed in the last few years to answer this question. 




FIG. 1: Pattern observed in erosion an experiment with com- 
mercial abrasive powder. Angles between 30° and 90° were 
found for different chevron alignments [3- 



For example, in Fig. [TJ we see a laboratory-scale experi- 
ment which reproduces a rich variety of natural patterns 
with few control parameters • These patterns are char- 
acterized by chevron alignments, what means that the 
paths created by the fluid present a characteristic angle 
that depends on the parameters varied in the experiment. 

In this paper, we investigate through numerical simu- 
lation a physical model that is designed to represent the 
generic situation of flowing water on a plane composed of 
credible sediment layer. This occurs naturally when the 
sea retreats from the shore or when a reservoir is drained. 
The flow of granular materials on planes is also of interest 
within the context of both industrial processing of pow- 
ders and geophysical instabilities such as landslides and 
avalanches. These flows have been found to be complex, 
exhibiting several different flow regimes as well as parti- 
cle segregation effects and instabilities [l^ . This paper is 
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FIG. 2: A typical initial configuration of the regularized ran- 
dom network with 64 grains and its corresponding triangula- 
tion. 



organized as follows. In Section II we describe the model 
used to simulate the erosion-deposition patterns. These 
patterns are presented and analysed in Section III, while 
the conclusions are left for Section IV. 



II. MODEL FORMULATION 

We model a system that takes into account the interac- 
tion of a granular medium with an incompressible Newto- 
nian fluid flowing through the corresponding pore space. 
The granular medium is initially considered as a, N x N 
regularized random network (RRN) on a plane where the 

Periodic boundaries 




FIG. 3: The velocity field is determined in the channel of 
length I and width w with one sphere of diameter d inside. 
The mean velocity of the fluid Vf is calculated in the rectan- 
gular cross-section. 



sites are the centers of mass of spherical grains with diam- 
eter d that are totally submerged in a fluid (e.g., water). 
We initialize this network as a regular square where the 
distance between the closest neighbors is Iq. The points 
(centers of mass of the grains) are then moved randomly 
along vectors with arbitrary direction and random mag- 
nitude that are smaller than the distance between the 
points. In this way, the points are distributed randomly, 
but with a characteristic distance. More precisely, in or- 
der to avoid the occurrence of overlapping grains, the 
maximum value adopted for the modulus of these dislo- 
cation vectors is {lo—d)/2. Although the lattice construc- 
tion is made in 2D, we are actually describing one layer 
of a three-dimensional system, since the centers of mass 
of the particles still lay on a 2D plane but the grains 
are considered to be spheres. This association will en- 
able us to generate a network of capillaries representing 
a complex geometry of the pore space. At this point, 
the entire system is triangulated considering each grain 
as a vertex of a Voronoi tesselation. In Fig. [2] we show a 
typical initial RRN conflguration and its corresponding 
triangulation 

Next, we assume that the local pore geometry between 
each two nearest neighbors grains i and j in this lattice 
can be modeled as a capillary channel of length I (the dis- 
tance between barycenters of the corresponding adjacent 
triangles), height d and width w equal to the distance 
between their centers. As shown in Fig. [3l if we consider 
periodic boundary conditions (PBC) in the x-direction, 
such a channel should be equivalent to a parallelepiped 
of fluid containing in its center a solid sphere of diameter 
d (grain). 

Here the Navier-Stokes and continuity equations in the 
three-dimensional channel are solved using the commer- 
cial CFD software FLUENT We consider no-slip 
boundary conditions at the bottom and assume that the 
top surface is shear-stress free. A pressure gradient Ap 
is imposed between the two ends of the channel in the 
z-direction. Its magnitude is sufficiently small to en- 
sure viscous flow conditions, i.e., a low Reynolds num- 
ber regime of flow. We adopt a non-structured tetra- 
hedral mesh to discretize the channel and an upwind 
finite-difference scheme is set to perform the numeri- 
cal simulations. The steady-state velocity and pres- 
sure fields are calculated for different channel geome- 
tries by systematically varying the ratio w/d in the range 
1.05 < w/d < 1.9. 

In Figs.|3^-c we show three contour plots of the velocity 
field computed for a channel with porosity w/d — 1.05 at 
heights y/d — A, 2 and 3/4, respectively. Considering the 
low Reynolds number conditions used in the simulations, 
the flow in the channel can be characterized in terms of 
a permeability index k through the relation. 



K Ap 
fi Az ' 



(1) 



where fi is the viscosity of the fluid and w/ is the average 
flow velocity. In practical terms, once the velocity and 
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FIG. 4: Contour plots of the velocity field obtained numer- 
ically for w/d — 1.05 at three different transverse planes of 
the channel, namely, (a) y/d = 1/4, (b) y/d = 1/2 and (c) 
y/d = 3/4. The gray shades ranging from dark to light cor- 
respond to low and high velocity magnitudes, respectively. 



pressure fields are obtained for a channel with a given 
value of w/d, we can compute the mean velocity value 
Vf through a cross-section orthogonal to the flow in the 
system. By repeating this procedure for different val- 
ues of the overall pressure drop Ap, we first confirm the 
validity of the linear relationship Vf oc Ap as expected 
from Eq. ([1]), so that the permeability k can be directly 
calculated from the slope of the corresponding straight 
line. As shown in Fig. [5l the dependence of k on the ra- 
tio w/d can be fully described as k/kq = f{w/d) where 
f{w/d) is a fourth degree polynomial of w/d and kq is 
the permeability for a channel with unitary aspect ratio, 
w/d= 1. 

Once the local geometry and permeability of all cap- 
illaries in the system is determined, we proceed by ap- 
plying a constant pressure drop between the inlet and 
outlet, i.e., the top and bottom of the entire pore net- 
work. Periodic boundary conditions are assumed in the 
lateral direction of the network, and the following local 
mass conservation equations are imposed at each of their 
N nodes to allow water flow throughout the entire pore 
space: 



(2) 



where the index j runs over all the neighbor nodes of node 
i, gij = wdn/l is the hydraulic conductance of the pore 
and Pi and pj are the pressures at nodes i and j, respec- 



FIG. 5: Permeability versus the ratio w/d for the local 
channel configuration shown in Fig. @. Here is the per- 
meability for a channel with ratio w/d = 1. The solid line is 
the best fit to the data of the fourth degree polynomial used 
to interpolate the local permeability of the channels in the 
pore network model. 



lively. Equation ^ corresponds to a set of A^^ coupled 
linear algebraic equations that are solved in terms of the 
nodal pressure field by means of a standard subroutine 
for sparse matrices. From the pressures in the nodes, the 
velocity magnitude of the fluid in each capillary can be 
computed. 

After the velocity field in the pore network is calcu- 
lated, we allow each grain to move in the system. By as- 
suming that there is no friction on the ground and drag is 
the only relevant force acting on the particles we obtain, 



dhc 



37r/xd^(v} 



v) 



(3) 



where is the fluid velocity at the channel i, v is the 
particle velocity, and the vectorial sum on the right is 
taken over all n channels surrounding the particle. By 
straightforward integration of the equation of motion ([31), 
the velocity and displacement caused by the fluid drag to 
each grain in the system during a time interval At can 
be written as, 



Ax 



At 
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(4) 



(5) 



where C = 18/i/pd^, p is the density of the grain, and Vg 
is its velocity in the previous time step. 

In our simulations, the center of each grain is displaced 
by Ax using a time step At that is sufficiently small to 
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numerically ensure that the pattern evolution remains 
invariant when compared with results performed using 
even smaller time steps. The distance between the top 
and bottom lines of the system is kept constant with its 
center of mass moving according to the mean displace- 
ment of all grains. What is crucial for the steadiness of 
the pattern formation is that one grain cannot overlap 
with another grain. After computing the movement of 
all grains at each time step, the pore space is then modi- 
fied, and we repeat the calculation of the velocity field to 
move the grains again, and so on. The simulation stops 
when the system reaches the steady-state, i.e., when the 
geometry of the aggregate remains unchanged with time. 



III. RESULTS 

We performed simulations on a lattice with 40 x 40 
grains with diameter d = 30/im and density p = 
2.75g/cm^. These particles are surrounded by water, i.e., 
a fluid of viscosity /i = 10~'^Pa.s which flows at small 
Reynolds conditions {Re <C 1) through a pore space of 
porosity that can be varied in the range 0.6 < (j) < 0.8. 
For each value of porosity we obtain results for five dif- 
ferent realizations of the initial random pore space. 

In Figs. [6^-c we show three final stable configurations 
of the model for different porosities values, 0=0.57, 0.71 
and 0.80, respectively. As can be observed, the steady- 
state patterns depend strongly on the porosity of the sys- 
tem. For sufficiently large values of the occurrence of 
particle clusters in the form of dendrites reflects the the 
strong coupling between fluid dynamics and grain move- 
ment, where the aligned preferential channels for flow 
leads to a high overall permeability of the porous system. 
For small porosities, however, no characteristic pattern is 
observed. This is to be expected since in compacted sys- 
tems the grains do not have much mobility while in loose 
systems the particles have the freedom to move in almost 
all directions. 

The tendency to form a dendritic pattern in which the 
particles align in preferential directions can be statisti- 
cally quantified if we determine for each pair of grains 
the angle a between the line connecting their centers of 
mass and the direction orthogonal to the flow (i.e., the 
x-direction shown in Fig. [S^). In Figs. [7^-c we show the 
histograms of these angles N{a) for 0=0.57, 0.71 and 
0.80. For a low porosity system, </> = 0.57, the results 
shown in Fig. [7^ indicate that a significant number of 
grains are aligned around a = 25° (and the symmet- 
ric direction of 155°), although the most frequent angle 
lies in the vicinity of 90°. This means that the parti- 
cles tend to be aligned in the vertical y-direction (i.e., 
the direction of the flux), what is expected as the system 
cannot change much from its initial configuration. In sys- 
tems with intermediate porosity values, 0.65 < 4> < 0.75, 
there is a substantial change in the preferred angle of 
alignment, which is around a = 60° (and the symmetric 
direction of 120°). This behaviour, as exemplified here 



in Fig. [Tjj, resembles the chevron alignment reported in 
Ref. where the results of the angle histograms are 
presented for a porosity = 0.71. In Fig. [71; we show the 
histogram N{a) for a large porosity system, 4> — 0.80, 
where no evident preferential direction can be observed, 
with the particles aligning themselves in angles between 
50° and 130°, with approximately the same probability. 

Finally, it is important to investigate the flow prop- 
erties of the porous system in terms of its macroscopic 
permeability Kj clS Oil function of porosity. In Fig. [8] we 
show that, while the permeability increases very slowly 
with porosity for diluted systems (i.e., for low values), 
it augments in a rather sharp way at high values of cj) to 
finally reach a divergent-like behavior in the vicinity of a 
porosity close to = 0.8. Interestingly, we find that this 
behavior can be well described by an inverse-logarithmic 
relation in the form k/kq = (a+6 In 0)~^ , as also depicted 
in Fig. m 



IV. CONCLUSIONS 

We developed a numerical model to describe the pro- 
cess of erosion-deposition caused by laminar flow with 
the drag force given by the Stokes law. The results we 
obtained show the formation of a typical erosion pattern 
characterized by chevron alignments very similar to the 
experimental ones presented in Ref. (a typical pattern 
is shown in Fig. [1]) . Through computational simulations 
performed with this model, we were able to find den- 
dritic patterns as well as to reproduce the preferential 
alignment of the chevron structures observed in real ex- 
periments. 

Our results indicate that these patterns depend sub- 
stantially on the porosity of the system. Previous stud- 
ies [U E, [l^ have shown that the size and shape of the 
particles influence dramatically the propagation of the 
fluid and the stress distribution in the system. Besides, 
it has been shown experimentally that the pattern ge- 
ometry must also depend on the flow properties through 
the porous medium |20l . [2l| , namely, on whether or not 
the inertial mechanisms of momentum transport play an 
important role on the dynamics of pattern formation. In 
the present study we considered only laminar flow. By 
changing the exponent of the velocity in the drag law, 
for example, one can reproduce aspects of turbulent flow 
to increase the complexity in the movement of the par- 
ticles. This could reveal a variety of new patterns. How 
the shape and the size distribution of the grains as well 
as the flow characteristics affect the patterns are natural 
questions that will be addressed in a future work. 

In recent works [ll|, [l^l, many patterns were ob- 
served in experiments involving avalanches, where one of 
the most important parameters is the depth of the sub- 
strate. Although the system studied here is related with 
erosion-sedimentation processes, this suggests that a va- 
riety of different patterns may be obtained just when one 
attempts to simulate them in three dimensions. A sim- 
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pic approximation to a three-dimensional system would 
be to consider that, depending on the flux, a particle 
would not stop as it reaches another particle, but could 
jump over it. With this possibility, the dynamics of the 
particles changes, as the velocity of them now depends 
also on the height of their centers of mass. 



and S. McNamara. This work was supported by the 
Brazilian agencies CNPq, CAPES and FUNCAP and 
Deutscher Akademischer Austauschdienst (DAAD). 
H. J. Herrmann thanks the Max-Planck prize. 



V. ACKNOWLEDGEMENTS 



We appreciate helpful interactions with 
A. M. C. Souza, A. A. P. Olarte, A. A. Moreira, 



[1] R. Behringer, H. Jaeger, and S. Nagel, Chaos 9, 509 
(1999). 

[2] H. M. Jaeger, S. R. Nagcl, and R. P. Behringer, Reviews 

of Modern Physics 68, 1259 (1996). 
[3] K. Kroy, G. Sauermann, and H. J. Herrmann, Phys. Rev. 

Lett. 88, 054301 (2002). 
[4] H. J. Herrmann, Physica A 263, 51 (1999). 
[5] D. Ertas and T. C. Halscy, Europh. Lett. 60, 931 (2002). 
[6] H. Seybold, J. S. Andrade, and H. J. Herrmann, PNAS 

104, 16804 (2007). 
[7] A. Daerr, P. Lee, J. Lanuza, and E. Clement, Phys. Rev. 

E 67, 065201 (R) (2003). 
[8] O. Pouliqucn, Physics of Fluids 11, 542 (1999). 
[9] L S. Aranson, F. Malloggi, and E. Clement, Phys. Rev. 

E 73, 050302(R) (2006). 
[10] O. Pouliquen and J. W. Vallance, Chaos 9, 621 (1999). 
[11] T. Borzsonyi, T. C. Halsey, and R. E. Ecke, Phys. Rev. 

Lett. 94, 208001 (2005). 
[12] Y. Forterre and O. PouUquen, J. Fluid Mech. 486, 21 

(2003). 

[13] O. Pouliquen, Phys. Rev. Lett. 93, 248001 (2004). 



[14] G. D. R. Midi, Euro Phys. J. E 14, 341 (2004). 

[15] O. Pouliquen, J. Dolour, and S. B. Savage, Letters to 
Nature 386, 816 (1997). 

[16] The FLUENT (trademark of FLUENT Inc.) fluid dy- 
namics analysis package has been used in this study; 
http: / /www. fluent. com. 

[17] J. Gong, D. HowcU, E. Longhi, R. P. Behringer, G. Rey- 
dcUct, L. Vancl, E. Clement, and S. Luding, Phys. Rev. 
Lett. 87, 035506 (2001). 

[18] G. Reydellet and E. Clement, Phys. Rev. Lett. 86, 3308 
(2001). 

[19] G. W. Baxter, R. P. Behringer, T. Fagert, and G. A. 

Johnson, Phys. Rev. Lett. 62, 2825 (1989). 
[20] J. S. Andrade, D. A. Street, T. Shinohara, Y. Shibusa, 

and Y. Aral, Phys. Rev. E 51, 5725 (1995). 
[21] J. S. Andrade, U. M. S. Costa, M. P. Almeida, H. A. 

Makse, and H. E. Stanley, Phys. Rev. Lett. 82, 5249 

(1999). 

[22] F. MaUoggi, J. Lanuza, B. Andreotti, and E. Clement, 
Europh. Lett. 75, 825 (2006). 



FIG. 6: In (a), (b) and (c) wc show the final configurations 
of the system for porosities 0=0.57, 0.71 and 0.80, respec- 
tively after 11000, 30000 and 140000 time steps. The grains 
arrange themselves in positions that give rise to characteristic 
patterns. 
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FIG. 7: Histograms of pairs of grains that arc touching and 
for which the line joining their centers of mass forms an angle 
a with the axis x. These results represent the average over 
five realizations of different porosities, namely, (a) <p = 0.57, 
(b) 0.71 and (c) 0.80. 



8 

80 \ 1 1 1 1 ' 1 ' — I ' 



60 




FIG. 8: Dependence of the macroscopic permeability k at 

the steady-state on the porosity (f> of the erosion-deposition 
system. The solid line is the best fit to the data of the function 
k/ko = {a + 61n^)~^, with parameters a = —0.21 and h = 
-1.01. 



